The present invention relates to processing of seismic data representative of subsurface features in the earth and, more particularly, to improved methods of and apparatus for processing seismic data using improved Radon transformations to remove unwanted noise from meaningful reflection signals.
Seismic surveys are one of the most important techniques for discovering the presence of oil and gas deposits. If the data are properly processed and interpreted, a seismic survey can give geologists a picture of subsurface geological features, so that they may better identify those features capable of holding oil and gas. Drilling is extremely expensive, and ever more so as easily tapped reservoirs are exhausted and new reservoirs are harder to reach. Having an accurate picture of an area's subsurface features can increase the odds of hitting an economically recoverable reserve and decrease the odds of wasting money and effort on a nonproductive well.
The principle behind seismology is deceptively simple. As seismic waves travel through the earth, portions of that energy are reflected back to the surface as the energy waves traverse different geological layers. Those seismic echoes or reflections give valuable information about the depth and arrangement of the formations, some of which hopefully contain oil or gas deposits.
A seismic survey is conducted by deploying an array of energy sources and an array of sensors or receivers in an area of interest. Typically, dynamite charges are used as sources for land surveys, and air guns are used for marine surveys. The sources are discharged in a predetermined sequence, sending seismic energy waves into the earth. The reflections from those energy waves or “signals” then are detected by the array of sensors. Each sensor records the amplitude of incoming signals over time at that particular location. Since the physical location of the sources and receivers is known, the time it takes for a reflection signal to travel from a source to a sensor is directly related to the depth of the formation that caused the reflection. Thus, the amplitude data from the array of sensors can be analyzed to determine the size and location of potential deposits.
This analysis typically starts by organizing the data from the array of sensors into common geometry gathers. That is, data from a number of sensors that share a common geometry are analyzed together. A gather will provide information about a particular spot or profile in the area being surveyed. Ultimately, the data will be organized into many different gathers and processed before the analysis is completed and the entire survey area mapped.
The types of gathers typically used include: common midpoint, where the sensors and their respective sources share a common midpoint; common source, where all the sensors share a common source; common offset, where all the sensors and their respective sources have the same separation or “offset”; and common receiver, where a number of sources share a common receiver. Common midpoint gathers are the most common gather today because they allow the measurement of a single point on a reflective subsurface feature from multiple source-receiver pairs, thus increasing the accuracy of the depth calculated for that feature.
The data in a gather are typically recorded or first assembled in the offset-time domain. That is, the amplitude data recorded at each of the receivers in the gather are assembled or displayed together as a function of offset, i.e., the distance of the receiver from a reference point, and as a function of time. The time required for a given signal to reach and be detected by successive receivers is a function of its velocity and the distance traveled. Those functions are referred to as kinematic travel time trajectories. Thus, at least in theory, when the gathered data are displayed in the offset-time domain, or “X-T” domain, the amplitude peaks corresponding to reflection signals detected at the gathered sensors should align into patterns that mirror the kinematic travel time trajectories. It is from those trajectories that one ultimately may determine an estimate of the depths at which formations exist.
A number of factors, however, make the practice of seismology and, especially, the interpretation of seismic data much more complicated than its basic principles. First, the reflected signals that indicate the presence of geological strata typically are mixed with a variety of noise.
The most meaningful signals are the so-called primary reflection signals, those signals that travel down to the reflective surface and then back up to a receiver. When a source is discharged, however, a portion of the signal travels directly to receivers without reflecting off of any subsurface features. In addition, a signal may bounce off of a subsurface feature, bounce off the surface, and then bounce off the same or another subsurface feature, one or more times, creating so-called multiple reflection signals. Other portions of the signal turn into noise as part of ground roll, refractions, and unresolvable scattered events. Some noise, both random and coherent, is generated by natural and man-made events outside the control of the survey.
All of this noise is occurring simultaneously with the reflection signals that indicate subsurface features. Thus, the noise and reflection signals tend to overlap when the survey data are displayed in X-T space. The overlap can mask primary reflection signals and make it difficult or impossible to identify patterns in the display upon which inferences about subsurface geology may be drawn. Accordingly, various mathematical methods have been developed to process seismic data in such a way that noise is separated from primary reflection signals.
Many such methods seek to achieve a separation of primary reflection signals and noise by transforming the data from the X-T domain to other domains. In other domains, such as the frequency-wavenumber (F-K) domain or the time-slowness (tau-P), there is less overlap between meaningful signal and noise data. Once the data are transformed, various mathematical filters are applied to the transformed data to separate noise from primary reflection signals, for example, by enhancing primary reflection signals or suppressing noise. The data then are inverse transformed back into the offset-time domain for further processing or interpretation.
For example, so-called Radon filters are commonly used to attenuate or remove multiple reflection signals. Such methods rely on Radon transformation equations to transform data from the offset-time (X-T) domain to the time-slowness (tau-P) where it can be filtered. More specifically, the X-T data are transformed along kinematic travel time trajectories having constant velocities and slownesses, where slowness p is defined as reciprocal velocity (or p=1/ν).
Radon filters have been developed for use in connection with common source, common receiver, common offset, and common midpoint gathers. They include those based on linear slant-stack, parabolic, and hyperbolic kinematic travel time trajectories. The general case forward transformation equation used in Radon filtration processes, R(p,τ)[d(x,t)], is set forth below:
            u      ⁡              (                  p          ,          τ                )              =                  ∫                  -          ∞                ∞            ⁢                        ⅆ          x                ⁢                              ∫                          -              ∞                        ∞                    ⁢                                    ⅆ              t                        ⁢                          ⅆ                              (                                  x                  ,                  t                                )                                      ⁢                          δ              ⁡                              [                                  f                  ⁡                                      (                                          t                      ,                      x                      ,                      τ                      ,                      p                                        )                                                  ]                                                              (          forward      ⁢                          ⁢      transformation        )  where                u(p,τ)=transform coefficient at slowness p and zero-offset time τ        d(x,t)=measured seismogram at offset x and two-way time t        p=slowness        t=two-way travel time        τ=two-way intercept time at p=0        x=offset        δ=Dirac delta function        ƒ(t,x,τ,p)=forward transform function        
The forward transform function for linear slant stack kinematic travel time trajectories is as follows:
            f      ⁡              (                  t          ,          x          ,          τ          ,          p                )              =          t      -      τ      -      px        where                                          δ            ⁡                          [                              f                ⁡                                  (                                      t                    ,                    x                    ,                    τ                    ,                    p                                    )                                            ]                                =                    ⁢                      δ            ⁡                          (                              t                -                τ                -                px                            )                                                                                =                        ⁢            1                    ,                                    when              ⁢                                                          ⁢              t                        =                          τ              +              px                                ,          and                                                          =                        ⁢            0                    ,                      elsewhere            .                              Thus, the forward linear slant stack Radon transformation equation becomes
      u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        x            ⁢              ⅆ                  (                      x            ,                          τ              +              px                                )                    
The forward transform function for parabolic kinematic trajectories is as follows:
            f      ⁡              (                  t          ,          x          ,          τ          ,          p                )              =          t      -      τ      -              px        2              where                                          δ            ⁡                          [                              f                ⁡                                  (                                      t                    ,                    x                    ,                    τ                    ,                    p                                    )                                            ]                                =                    ⁢                      δ            ⁡                          (                              t                -                τ                -                                  px                  2                                            )                                                                                =                        ⁢            1                    ,                                    when              ⁢                                                          ⁢              t                        =                          τ              +                              px                2                                              ,          and                                                          =                        ⁢            0                    ,                      elsewhere            .                              Thus, the forward parabolic Radon transformation equation becomes
      u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        x            ⁢              ⅆ                  (                      x            ,                          τ              +                              px                2                                              )                    
The forward transform function for hyperbolic kinematic travel time trajectories is as follows:
            f      ⁡              (                  t          ,          x          ,          τ          ,          p                )              =          t      -                                    τ            2                    +                                    p              2                        ⁢                          x              2                                            where                                          δ            ⁡                          [                              f                ⁡                                  (                                      t                    ,                    x                    ,                    τ                    ,                    p                                    )                                            ]                                =                    ⁢                      δ            ⁡                          (                              t                -                                                                            τ                      2                                        +                                                                  p                        2                                            ⁢                                              x                        2                                                                                                        )                                                                                =                        ⁢            1                    ,                                    when              ⁢                                                          ⁢              t                        =                                                            τ                  2                                +                                                      p                    2                                    ⁢                                      x                    2                                                                                ,          and                                                          =                        ⁢            0                    ,                      elsewhere            .                              Thus, the forward hyperbolic Radon transformation equation becomes
      u    ⁡          (              p        ,        τ            )        =            ∫              -        ∞            ∞        ⁢                  ⅆ        x            ⁢              ⅆ                  (                      x            ,                                                            τ                  2                                +                                                      p                    2                                    ⁢                                      x                    2                                                                                )                    
A general case inverse Radon transformation equation is set forth below:
            d      ⁡              (                  x          ,          t                )              =                  ∫                  -          ∞                ∞            ⁢                        ⅆ          p                ⁢                              ∫                          -              ∞                        ∞                    ⁢                                    ⅆ                              τρ                ⁡                                  (                  τ                  )                                                      *                          u              ⁡                              (                                  p                  ,                  τ                                )                                      ⁢                          δ              ⁡                              [                                  g                  ⁡                                      (                                          t                      ,                      x                      ,                      τ                      ,                      p                                        )                                                  ]                                                              (          inverse      ⁢                          ⁢      transformation        )  where                g(t,x,τ,p)=inverse transform function, and        ρ(τ)*=convolution of rho filter.        
The inverse transform function for linear slant stack kinematic trajectories is as follows:g(t,x,τ,p)=τ−t+px Thus, when τ=t−px, the inverse linear slant stack Radon transformation equation becomes
      d    ⁡          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢                  ⁢                  ⅆ        p            ⁢                          ⁢              ρ        ⁡                  (          τ          )                    *              u        ⁡                  (                      p            ,                          t              -              px                                )                    
The inverse transform function for parabolic trajectories is as follows:g(t,x,τ,p)=τ−t+px2 Thus, when τ=t−px2, the inverse linear slant stack Radon transformation equation becomes
      d    ⁡          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢                  ⁢                  ⅆ        p            ⁢                          ⁢              ρ        ⁡                  (          τ          )                    *              u        ⁡                  (                      p            ,                          t              -                              px                2                                              )                    
The inverse transform function for hyperbolic trajectories is as follows:g(t,x,τ,p)=τ−√{square root over (t2−p2x2)}Thus, when τ=√{square root over (t2−p2x2)}, the inverse hyperbolic Radon transformation equation becomes
      d    ⁡          (              x        ,        t            )        =            ∫              -        ∞            ∞        ⁢                  ⁢                  ⅆ        p            ⁢                          ⁢              ρ        ⁡                  (          τ          )                    *              u        ⁡                  (                      p            ,                                                            t                  2                                -                                                      p                    2                                    ⁢                                      x                    2                                                                                )                    
The choice of which form of Radon transformation preferably is guided by the travel time trajectory at which signals of interest in the data are recorded. Common midpoint gathers, because they offer greater accuracy by measuring a single point from multiple source-receiver pairs, are preferred over other types of gathers. Primary reflection signals in a common midpoint gather generally will have hyperbolic kinematic trajectories. Thus, it would be preferable to use hyperbolic Radon transforms.
To date, however, Radon transformations based on linear slant stack and parabolic trajectories have been more commonly used. The transform function for hyperbolic trajectories contains a square root whereas those for linear slant stack and parabolic transform functions do not. Thus, the computational intensity of hyperbolic Radon transformations has in large part discouraged their use in seismic processing.
It has not been practical to accommodate the added complexity of hyperbolic Radon transformations because the computational requirements of conventional processes are already substantial. For example, such prior art Radon methods typically first process the data to compensate for the increase in travel time as sensors are further removed from the source. This step is referred to as normal moveout or “NMO” correction. It is designed to eliminate the differences in time that exist between the primary reflection signals recorded at close-in receivers, i.e., at near offsets, and those recorded at remote receivers, i.e., at far offsets.
NMO correction involves a least-mean-squares (“LMS”) energy minimization calculation. Forward and inverse Radon transforms also are not exact inverses of each other. Accordingly, an additional LMS calculation, conjugate gradient and/or sparsity constraints are often used in the transformation. Those calculations in general and, especially LMS analyses, are complex and require substantial computing time, even on advanced computer systems.
Moreover, a typical seismic survey may involve hundreds of sources and receivers, thus generating tremendous quantities of raw data. The data may be subjected to thousands of different data gathers. Each gather is subjected not only to filtering processes as described above, but in all likelihood to many other enhancement processes as well. For example, data are typically processed to compensate for the diminution in amplitude as a signal travels through the earth (“amplitude balancing”). Then, once the individual gathers have been filtered, they must be “stacked”, or compiled with other gathers, and subjected to further processing in order to make inferences about subsurface formations. Seismic processing by prior art Radon methods, therefore, requires substantial and costly computational resources, and there is a continuing need for more efficient and economical processing methods.
For example, French Pat. 2,801,387 of B. Fernand (“Fernand”) relates generally to processing of seismic data and specifically to methods that utilize a determination of residual dynamic corrections (i.e., residual “normal move out” or “NMO” corrections). It generally discloses a method wherein amplitude data in the X-T domain are transformed into the tau-P domain using a parabolic Radon transformation. A filter is defined and applied in the tau-P domain which attenuates all data, including the transformed primary signals, below a defined slowness value pmute. The data that remain after attenuating the primary reflection signals contain a substantial portion of the signals corresponding to multiple reflections. The unmuted data are then transformed back into offset-time space and are subtracted from the original data in the gather. The subtraction process removes the multiple reflection signals from the data gather, leaving the primary reflection signals more readily apparent and easier to interpret. It will be appreciated, however, that such prior art Radon filters are computationally intense, in large part because they rely on NMO correction, and rely on less accurate parabolic Radon transformations instead of hyperbolic transformations.
The velocity at which primary reflection signals are traveling, what is referred to as the stacking velocity, is used in various analytical methods that are applied to seismic data. For example, it is used in determining the depth and lithology or sediment type of geological formations in the survey area. It also is used various seismic attribute analyses. A stacking velocity function may be determined by what is referred to as a semblance analysis. Which such methods have provided reasonable estimations of the stacking velocity function for a data set, given the many applications based thereon, it is desirable to define the stacking velocity function as precisely as possible.
An object of this invention, therefore, is to provide a method for processing seismic data that more effectively removes unwanted noise from meaningful reflection signals.
It also is an object to provide such methods based on hyperbolic Radon transformations.
Another object of this invention is to provide methods for removal of noise from seismic data that are comparatively simple and require relatively less computing time and resources.
Yet another object is to provide methods for more accurately defining the stacking velocity function for a data set.
It is a further object of this invention to provide such methods wherein all of the above-mentioned advantages are realized.
Those and other objects and advantages of the invention will be apparent to those skilled in the art upon reading the following detailed description and upon reference to the drawings.